knitr::opts_chunk$set(echo = TRUE)

rm(list = ls())

#cmdstanr::set_cmdstan_path(path = "C:/Users/kueng/.cmdstan/cmdstan-2.35.0")
#cmdstanr::set_cmdstan_path(path = "C:/Users/pascku/.cmdstan/cmdstan-2.36.0")

library(tidyverse)
library(R.utils)
library(wbCorr)
library(readxl)
library(kableExtra)
library(brms)
library(bayesplot)
library(see)
library(beepr)
library(DHARMa)
library(digest)



source(file.path('Functions', 'ReportModels.R'))
source(file.path('Functions', 'PrettyTables.R'))
source(file.path('Functions', 'ReportMeasures.R'))
source(file.path('Functions', 'PrepareData.R'))

report_function_hash <- digest::digest(summarize_brms)
system("shutdown /a")
## [1] 1116
# Set options for analysis
use_mi = FALSE
shutdown = FALSE
report_ordinal = FALSE
do_priorsense = FALSE
get_bayesfactor = TRUE
check_models = TRUE 

if (get_bayesfactor) {
  stats_to_report <- c('CI', 'SE', 'pd', 'ROPE', 'BF', 'Rhat', 'ESS')
} else {
  stats_to_report <- c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS')
}

options(
  dplyr.print_max = 100, 
  brms.backend = 'cmdstan',
  brms.file_refit = ifelse(use_mi, 'never', 'on_change'),
  brms.file_refit = 'on_change',
  #brms.file_refit = 'always',
  error = function() {
    beepr::beep(sound = 5)
    if (shutdown) {
      system("shutdown /s /t 180")
      quit(save = "no", status = 1)
    }
  }
  , es.use_symbols = TRUE
)


####################### Model parameters #######################

iterations = 12000 # 12'000 per chain to achieve 40'000
warmup = 2000 # 2000

# NO AR!!!
#corstr = 'ar'
#corstr = 'cosy_couple'
#corstr = 'cosy_couple:user'


################################################################

suffix = paste0('_final_with_plan_', as.character(iterations))
df <- openxlsx::read.xlsx(file.path('long.xlsx'))
df_original <- df

df_double <- prepare_data(df, recode_pushing = TRUE, use_mi = use_mi)[[1]]

Constructing scales Re-coding pusing reshaping data (4field) centering data within and between

summary(df_double$pushing)

Min. 1st Qu. Median Mean 3rd Qu. Max. NA’s 0.0000 0.0000 0.0000 0.1649 0.0000 5.0000 275

Modelling

# For indistinguishable Dyads
model_rows_fixed <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'plan_selfPlan',
    'plan_partnerPlan',
    'weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'pressure_self_cb',
    'pressure_partner_cb',
    'pushing_self_cb',
    'pushing_partner_cb',
    'weartime_self_cb'
  )


model_rows_fixed_ordinal <- c(
  model_rows_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rows_fixed[2:length(model_rows_fixed)]
)

model_rows_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rows_random_ordinal <- c(model_rows_random,'disc')
# For indistinguishable Dyads
model_rownames_fixed <- c(
    "Intercept", 
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Daily individual's experienced persuasion",  
    "Daily partner's experienced persuasion", 
    "Daily individual's experienced pressure", 
    "Daily partner's experienced pressure", 
    "Daily individual's experienced pushing", 
    "Daily partner's experienced pushing", 
    "Day", 
    "Own Actionplan",
    'Partner Actionplan',
    "Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Mean individual's experienced persuasion", 
    "Mean partner's experienced persuasion", 
    "Mean individual's experienced pressure", 
    "Mean partner's experienced pressure", 
    "Mean individual's experienced pushing", 
    "Mean partner's experienced pushing", 
    "Mean weartime"
  )


model_rownames_fixed_ordinal <- c(
  model_rownames_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed[2:length(model_rownames_fixed)]
)

model_rownames_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  "sd(Daily individual's experienced persuasion)", 
  "sd(Daily partner's experienced persuasion)", # OR partner received
  "sd(Daily individual's experienced pressure)", 
  "sd(Daily partner's experienced pressure)", 
  "sd(Daily individual's experienced pushing)", 
  "sd(Daily partner's experienced pushing)", 
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rownames_random_ordinal <- c(model_rownames_random,'disc')
rows_to_pack <- list(
  "Within-Person Effects" = c(2,11),
  "Between-Person Effects" = c(12,18),
  "Random Effects" = c(19, 25), 
  "Additional Parameters" = c(26,26)
  )


rows_to_pack_ordinal <- list(
  "Within-Person Effects" = c(2+5,11+5),
  "Between-Person Effects" = c(12+5,18+5),
  "Random Effects" = c(19+5, 25+5), 
  "Additional Parameters" = c(26+5,26+6)
  )

Self-Reported MVPA

range(df_double$pa_sub, na.rm = T) 
## [1]   0 720
hist(df_double$pa_sub, breaks = 40) 

hist(log(df_double$pa_sub+00000000001), breaks = 40)

Hurdle Lognormal Model

formula <- bf(
  pa_sub ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    plan_self + plan_partner +
    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | dd | coupleID),
  
  hu = ~ persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    plan_self + plan_partner +
    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | dd | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
) 

prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
  , brms::set_prior("normal(0, 50)", class = "Intercept") # for non-zero PA
  , brms::set_prior("normal(0.5, 2.5)", class = "Intercept", dpar = 'hu') # hurdle part
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = hurdle_lognormal()
#)

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_sub <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::hurdle_lognormal(), 
  #family = brms::hurdle_negbinomial(), 
  #family = brms::hurdle_poisson(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 42,
  file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
pa_sub_digest <- digest::digest(pa_sub)
if (check_models) {
  check_brms(pa_sub, log_pp_check = TRUE)
  DHARMa.check_brms.all(pa_sub, integer = TRUE, outliers_type = 'bootstrap')
}
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                   Term  VIF   VIF 95% CI Increased SE Tolerance
##     persuasion_self_cw 1.59 [1.54, 1.65]         1.26      0.63
##  persuasion_partner_cw 1.08 [1.05, 1.11]         1.04      0.93
##       pressure_self_cw 1.07 [1.04, 1.10]         1.03      0.94
##    pressure_partner_cw 1.05 [1.03, 1.09]         1.03      0.95
##        pushing_self_cw 1.04 [1.02, 1.08]         1.02      0.96
##     pushing_partner_cw 1.09 [1.06, 1.12]         1.04      0.92
##     persuasion_self_cb 1.08 [1.06, 1.12]         1.04      0.92
##  persuasion_partner_cb 3.68 [3.52, 3.86]         1.92      0.27
##       pressure_self_cb 3.68 [3.52, 3.86]         1.92      0.27
##    pressure_partner_cb 2.23 [2.14, 2.32]         1.49      0.45
##        pushing_self_cb 2.19 [2.10, 2.28]         1.48      0.46
##     pushing_partner_cb 3.14 [3.01, 3.29]         1.77      0.32
##              plan_self 3.09 [2.96, 3.23]         1.76      0.32
##           plan_partner 1.40 [1.36, 1.45]         1.18      0.71
##                    day 1.37 [1.33, 1.42]         1.17      0.73
##  Tolerance 95% CI
##      [0.60, 0.65]
##      [0.90, 0.95]
##      [0.91, 0.96]
##      [0.92, 0.97]
##      [0.92, 0.98]
##      [0.89, 0.94]
##      [0.89, 0.95]
##      [0.26, 0.28]
##      [0.26, 0.28]
##      [0.43, 0.47]
##      [0.44, 0.48]
##      [0.30, 0.33]
##      [0.31, 0.34]
##      [0.69, 0.74]
##      [0.70, 0.75]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        cauchy        100%
## 
## Predicted Distribution of Response
## 
##                Distribution Probability
##  neg. binomial (zero-infl.)         78%
##               beta-binomial         16%
##                   lognormal          3%
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 1 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 9, observations = 3684, p-value = 0.34
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.0005428882 0.0029858849
## sample estimates:
## outlier frequency (expected: 0.0016042345276873 ) 
##                                       0.002442997
if (do_priorsense) {
  priorsense_vars <- c(
      'Intercept',
      'b_persuasion_self_cw',
      'b_persuasion_partner_cw',
      'b_pressure_self_cw',
      'b_pressure_partner_cw',
      'b_pushing_self_cw',
      'b_pushing_partner_cw'
  )
  
  hurdle_priorsense_vars <- c(
    'Intercept_hu',
    'b_hu_persuasion_self_cw',
    'b_hu_persuasion_partner_cw',
    'b_hu_pressure_self_cw',
    'b_hu_pressure_partner_cw',
    'b_hu_pushing_self_cw',
    'b_hu_pushing_partner_cw'
  )
  
  gc()
  priorsense::powerscale_sensitivity(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
  priorsense::powerscale_plot_dens(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
  priorsense::powerscale_plot_ecdf(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
  priorsense::powerscale_plot_quantities(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
}
# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub$data$pa_sub[pa_sub$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_sub <- summarize_brms(
  pa_sub, 
  stats_to_report = stats_to_report,
  rope_range = rope_range_continuous,
  hu_rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) 
## Sampling priors, please wait...
## Warning in summarize_brms(pa_sub, stats_to_report = stats_to_report, rope_range
## = rope_range_continuous, : Coefficients were exponentiated. Double check if
## this was intended.
# Print the updated dataframe
summary_pa_sub %>%
  print_df(rows_to_pack = rows_to_pack)
exp(Est.)_hu SE_hu 95% CI_hu pd_hu ROPE_hu inside ROPE_hu BF_hu BF_Evidence_hu Rhat_hu Bulk_ESS_hu Tail_ESS_hu exp(Est.)_nonzero SE_nonzero 95% CI_nonzero pd_nonzero ROPE_nonzero inside ROPE_nonzero BF_nonzero BF_Evidence_nonzero Rhat_nonzero Bulk_ESS_nonzero Tail_ESS_nonzero
Intercept 0.21*** 0.04 [ 0.15, 0.30] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 46489 34693 35.17*** 2.95 [29.76, 41.52] 1.000 [0.92, 1.08] 0.000 >100 Overwhelming Evidence 1.000 17579 26387
Within-Person Effects
Daily individual’s experienced persuasion 1.58*** 0.11 [ 1.37, 1.84] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 50906 32990 1.03 0.03 [ 0.98, 1.09] 0.869 [0.92, 1.08] 0.961 0.018 Very Strong Evidence for Null 1.000 35141 34015
Daily partner’s experienced persuasion 1.34*** 0.09 [ 1.18, 1.53] 1.000 [0.84, 1.20] 0.036 >100 Overwhelming Evidence 1.000 56832 32154 1.03 0.02 [ 0.99, 1.08] 0.910 [0.92, 1.08] 0.977 0.020 Very Strong Evidence for Null 1.000 43176 35087
Daily individual’s experienced pressure 0.98 0.16 [ 0.69, 1.36] 0.556 [0.84, 1.20] 0.727 0.080 Strong Evidence for Null 1.000 55722 31290 0.90* 0.05 [ 0.80, 0.99] 0.980 [0.92, 1.08] 0.277 0.070 Strong Evidence for Null 1.000 45139 39256
Daily partner’s experienced pressure 1.55* 0.34 [ 1.04, 2.67] 0.983 [0.84, 1.20] 0.103 1.065 Weak Evidence 1.000 43532 28265 0.95 0.04 [ 0.86, 1.03] 0.894 [0.92, 1.08] 0.688 0.014 Very Strong Evidence for Null 1.000 42665 36112
Daily individual’s experienced pushing 1.01 0.14 [ 0.76, 1.37] 0.518 [0.84, 1.20] 0.791 0.069 Strong Evidence for Null 1.000 47778 32843 0.99 0.03 [ 0.93, 1.06] 0.609 [0.92, 1.08] 0.971 0.008 Very Strong Evidence for Null 1.000 49766 34692
Daily partner’s experienced pushing 1.34** 0.16 [ 1.07, 1.72] 0.995 [0.84, 1.20] 0.165 1.658 Weak Evidence 1.000 56748 32068 0.96 0.03 [ 0.90, 1.02] 0.892 [0.92, 1.08] 0.891 0.015 Very Strong Evidence for Null 1.000 52279 34615
Day 0.81 0.12 [ 0.60, 1.09] 0.914 [0.84, 1.20] 0.423 0.194 Moderate Evidence for Null 1.000 91018 29772 0.98 0.06 [ 0.87, 1.11] 0.599 [0.92, 1.08] 0.779 0.007 Very Strong Evidence for Null 1.000 80027 30660
Own Actionplan 10.50*** 1.12 [ 8.57, 12.90] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 82098 33249 1.34*** 0.06 [ 1.21, 1.47] 1.000 [0.92, 1.08] 0.000 >100 Overwhelming Evidence 1.000 80029 30755
Partner Actionplan 1.17 0.12 [ 0.95, 1.43] 0.933 [0.84, 1.20] 0.602 0.154 Moderate Evidence for Null 1.000 80829 32856 1.07 0.05 [ 0.98, 1.17] 0.937 [0.92, 1.08] 0.599 0.024 Very Strong Evidence for Null 1.000 86439 28993
Daily weartime NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion 1.41 0.54 [ 0.65, 3.02] 0.815 [0.84, 1.20] 0.242 0.293 Moderate Evidence for Null 1.000 29036 30407 0.96 0.15 [ 0.70, 1.30] 0.612 [0.92, 1.08] 0.373 0.015 Very Strong Evidence for Null 1.000 26670 29681
Mean partner’s experienced persuasion 1.32 0.51 [ 0.61, 2.80] 0.764 [0.84, 1.20] 0.282 0.250 Moderate Evidence for Null 1.000 28814 30891 0.90 0.14 [ 0.66, 1.22] 0.753 [0.92, 1.08] 0.310 0.019 Very Strong Evidence for Null 1.000 26031 30027
Mean individual’s experienced pressure 0.17*** 0.10 [ 0.05, 0.49] 1.000 [0.84, 1.20] 0.001 67.790 Very Strong Evidence 1.000 42758 33569 0.76 0.28 [ 0.36, 1.63] 0.769 [0.92, 1.08] 0.126 0.030 Strong Evidence for Null 1.000 10949 18226
Mean partner’s experienced pressure 0.25** 0.14 [ 0.08, 0.72] 0.996 [0.84, 1.20] 0.011 7.904 Moderate Evidence 1.000 42429 34022 0.56 0.21 [ 0.26, 1.20] 0.936 [0.92, 1.08] 0.050 0.076 Strong Evidence for Null 1.000 10613 17470
Mean individual’s experienced pushing 1.42 0.83 [ 0.45, 4.62] 0.727 [0.84, 1.20] 0.201 0.350 Weak Evidence for Null 1.000 33685 32275 1.51 0.41 [ 0.87, 2.62] 0.932 [0.92, 1.08] 0.073 0.048 Strong Evidence for Null 1.000 15413 24107
Mean partner’s experienced pushing 2.58 1.53 [ 0.81, 8.39] 0.945 [0.84, 1.20] 0.070 1.066 Weak Evidence 1.000 36276 33152 1.75* 0.49 [ 1.00, 3.04] 0.976 [0.92, 1.08] 0.031 0.118 Moderate Evidence for Null 1.000 14406 23651
Mean weartime NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.81 0.11 [0.62, 1.08] NA NA NA NA NA 1.000 28946 31750 0.30 0.04 [0.23, 0.40] NA NA NA NA NA 1.000 18892 27189
sd(Daily individual’s experienced persuasion) 0.22 0.09 [0.04, 0.41] NA NA NA NA NA 1.000 16555 13405 0.11 0.02 [0.07, 0.16] NA NA NA NA NA 1.000 29886 31572
sd(Daily partner’s experienced persuasion) 0.15 0.10 [0.01, 0.35] NA NA NA NA NA 1.000 15649 20545 0.08 0.02 [0.04, 0.13] NA NA NA NA NA 1.000 26632 24369
sd(Daily individual’s experienced pressure) 0.20 0.18 [0.01, 0.77] NA NA NA NA NA 1.000 25019 26499 0.07 0.06 [0.00, 0.23] NA NA NA NA NA 1.000 23208 25726
sd(Daily partner’s experienced pressure) 0.36 0.30 [0.02, 1.21] NA NA NA NA NA 1.000 18888 25903 0.06 0.05 [0.00, 0.18] NA NA NA NA NA 1.000 26047 22976
sd(Daily individual’s experienced pushing) 0.51 0.17 [0.20, 0.93] NA NA NA NA NA 1.000 22851 19964 0.09 0.04 [0.02, 0.17] NA NA NA NA NA 1.000 17279 15406
sd(Daily partner’s experienced pushing) 0.21 0.15 [0.01, 0.57] NA NA NA NA NA 1.000 20680 24410 0.08 0.04 [0.01, 0.16] NA NA NA NA NA 1.000 17582 15911
Additional Parameters
sigma NA NA NA NA NA NA NA NA NA NA NA 0.68 0.01 [0.65, 0.70] NA NA NA NA NA 1.000 69410 29347 
# Plot continuous part of model

variable <- c(
  '(Intercept)',
  'b_persuasion_self_cw',
  'b_persuasion_partner_cw',
  'b_pressure_self_cw',
  'b_pressure_partner_cw',
  'b_pushing_self_cw',
  'b_pushing_partner_cw'
)


plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(
    pa_sub, 
    parameter = variable, 
    range = rope_range_continuous,
    verbose = F,
    ci = 1
  )
) + theme_bw()

# Hurdle part of the model
variable <- c(
  'b_hu_persuasion_self_cw',
  'b_hu_persuasion_partner_cw',
  'b_hu_pressure_self_cw',
  'b_hu_pressure_partner_cw',
  'b_hu_pushing_self_cw',
  'b_hu_pushing_partner_cw'
)

plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length

# The rope range for the bernoulli part of the model is -0.18, 0.18
plot(
  bayestestR::rope(pa_sub, parameter = variable, range = c(-0.18, 0.18), ci = 1),
  verbose = FALSE
) + theme_bw()
## Possible multicollinearity between b_hu_persuasion_partner_cb and
##   b_hu_persuasion_self_cb (r = 0.77), b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.76), b_pressure_partner_cb and
##   b_pressure_self_cb (r = 0.82). This might lead to inappropriate results.
##   See 'Details' in '?rope'.

Hurdle part of the model on the left, non-zero part towards the right side of the table

conds_eff <- conditional_spaghetti(
  pa_sub, 
  effects = c(
    'persuasion_self_cw',
    'persuasion_partner_cw',
    'pressure_self_cw',
    'pressure_partner_cw',
    'pushing_self_cw',
    'pushing_partner_cw'
  ),
  x_label = c(
    'Received Persuasion',
    'Exerted Persuasion',
    'Received Pressure',
    'Exerted Pressure',
    'Received Plan-Related Pushing',
    'Exerted Plan-Related Pushing'
  ),
  group_var = 'coupleID',
  plot_full_range = TRUE,
  y_limits = c(0, 100),
  y_label = "Same-Day MVPA",
  y_labels = c('Probability of Being Active', 'Minutes of MVPA When Active', 'Overall Expected Minutes of MVPA'),
  , filter_quantiles = .9995
  , font_family = 'Candara'
)
## This is posterior version 1.6.0
## 
## Attaching package: 'posterior'
## The following object is masked from 'package:bayesplot':
## 
##     rhat
## The following objects are masked from 'package:stats':
## 
##     mad, sd, var
## The following objects are masked from 'package:base':
## 
##     %in%, match
## Registering fonts with R
## Scanning ttf files in C:\Windows/Fonts, C:\Users\pascku\AppData\Local/Microsoft/Windows/Fonts ...
## Extracting .afm files from .ttf files...
## C:\Windows\Fonts\Candara.ttf : Candara already registered in fonts database. Skipping.
## C:\Windows\Fonts\Candarab.ttf : Candara-Bold already registered in fonts database. Skipping.
## C:\Windows\Fonts\Candarai.ttf : Candara-Italic already registered in fonts database. Skipping.
## C:\Windows\Fonts\Candaral.ttf : Candara-Light already registered in fonts database. Skipping.
## C:\Windows\Fonts\Candarali.ttf : Candara-LightItalic already registered in fonts database. Skipping.
## C:\Windows\Fonts\Candaraz.ttf : Candara-BoldItalic already registered in fonts database. Skipping.
## Found FontName for 0 fonts.
## Scanning afm files in C:/Users/pascku/AppData/Local/R/cache/R/renv/cache/v5/windows/R-4.4/x86_64-w64-mingw32/extrafontdb/1.0/a861555ddec7451c653b40e713166c6f/extrafontdb/metrics
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
print(conds_eff)

$persuasion_self_cw

## Warning: Removed 63 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 16 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.00845

$persuasion_partner_cw

## Warning: Removed 37 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.00661

$pressure_self_cw

## Picking joint bandwidth of 0.012

$pressure_partner_cw

## Picking joint bandwidth of 0.0222

$pushing_self_cw

## Picking joint bandwidth of 0.0107

$pushing_partner_cw

## Picking joint bandwidth of 0.011

Note. This graphic illustrates the relationship between social control and moderate to vigorous physical activity (MVPA) using a Bayesian Hurdle-Lognormal Multilevel Model. The predictor is centered within individuals to examine how deviations from their average social control relate to same-day MVPA. Shaded areas indicate credible intervals, thick lines show fixed effects, and thin lines represent random effects, highlighting variability across couples. The plots display the probability of being active, expected minutes of MVPA when active, and combined predicted MVPA. The bottom density plot visualizes the posterior distributions of slope estimates, transformed to represent multiplicative changes in odds ratios (hurdle component) or expected values. Medians and 95% credible intervals (2.5th and 97.5th percentiles) are shown. Effects are significant, when the 95% credible interval does not overlap 1.

x_label = c(
    'Received Persuasion',
    'Exerted Persuasion',
    'Received Pressure',
    'Exerted Pressure',
    'Received Plan-Related Pushing',
    'Exerted Plan-Related Pushing'
  )

home_dir <- getwd()
output_dir <- file.path(home_dir, 'Output', 'Plots')

for (i in 1:length(conds_eff)) {
  effname <- names(conds_eff)[i]
  eff_plot <- conds_eff[[i]]
  x_label_i <- x_label[[i]]
  
  rmarkdown::render(
    file.path(output_dir, 'BeautifulPlotWithNote.Rmd'), 
    output_file = file.path(output_dir, paste0('Graphic_', effname, '.pdf')),
    params = list(
      home_dir = home_dir,
      output_dir = output_dir,
      p_i = eff_plot,
      p_name = effname,
      x_label = x_label_i
      ),
    envir = new.env(),
    quiet = TRUE
  )
}

print('done')

Comparing effect size of pressure and pushing

hypothesis(pa_sub, "pressure_self_cw < pushing_self_cw")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (pressure_self_cw... < 0     -0.1      0.07    -0.21     0.01      13.68
##   Post.Prob Star
## 1      0.93     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Device Based MVPA

range(df_double$pa_obj, na.rm = T) 
## [1]   5.75 971.25
hist(df_double$pa_obj, breaks = 50)

df_double$pa_obj_log <- log(df_double$pa_obj)

hist(df_double$pa_obj_log, breaks = 50)

Lognormal Model

formula <- bf(
  pa_obj ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    plan_self + plan_partner +
    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 50)", class = "Intercept") 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = lognormal()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_obj_log <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = lognormal(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
pa_obj_log_digest <- digest::digest(pa_obj_log)
if (check_models) {
  check_brms(pa_obj_log, log_pp_check = TRUE)
  DHARMa.check_brms.all(pa_obj_log, integer = TRUE, outliers_type = 'bootstrap')
}
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                   Term  VIF    VIF 95% CI Increased SE Tolerance
##     persuasion_self_cw 1.11 [1.08,  1.15]         1.05      0.90
##  persuasion_partner_cw 1.13 [1.09,  1.17]         1.06      0.89
##       pressure_self_cw 1.07 [1.04,  1.12]         1.03      0.93
##    pressure_partner_cw 1.10 [1.07,  1.15]         1.05      0.91
##        pushing_self_cw 1.13 [1.09,  1.17]         1.06      0.89
##     pushing_partner_cw 1.18 [1.15,  1.23]         1.09      0.85
##       pressure_self_cb 3.96 [3.75,  4.18]         1.99      0.25
##    pressure_partner_cb 4.33 [4.10,  4.57]         2.08      0.23
##              plan_self 1.29 [1.24,  1.34]         1.14      0.78
##           plan_partner 1.30 [1.25,  1.35]         1.14      0.77
##                    day 1.03 [1.01,  1.09]         1.02      0.97
##       weartime_self_cw 1.04 [1.02,  1.09]         1.02      0.96
##       weartime_self_cb 1.25 [1.20,  1.30]         1.12      0.80
##  Tolerance 95% CI
##      [0.87, 0.93]
##      [0.85, 0.91]
##      [0.90, 0.96]
##      [0.87, 0.93]
##      [0.85, 0.92]
##      [0.81, 0.87]
##      [0.24, 0.27]
##      [0.22, 0.24]
##      [0.74, 0.80]
##      [0.74, 0.80]
##      [0.92, 0.99]
##      [0.92, 0.98]
##      [0.77, 0.83]
## 
## Moderate Correlation
## 
##                   Term  VIF    VIF 95% CI Increased SE Tolerance
##     persuasion_self_cb 9.09 [8.58,  9.64]         3.02      0.11
##  persuasion_partner_cb 9.90 [9.34, 10.50]         3.15      0.10
##        pushing_self_cb 6.40 [6.05,  6.78]         2.53      0.16
##     pushing_partner_cb 6.22 [5.88,  6.59]         2.49      0.16
##  Tolerance 95% CI
##      [0.10, 0.12]
##      [0.10, 0.11]
##      [0.15, 0.17]
##      [0.15, 0.17]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        cauchy         88%
##     lognormal          9%
##       weibull          3%
## 
## Predicted Distribution of Response
## 
##                Distribution Probability
##                   lognormal         72%
##                     tweedie         16%
##  neg. binomial (zero-infl.)          9%
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 27, observations = 3292, p-value < 2.2e-16
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.001663123 0.005467801
## sample estimates:
## outlier frequency (expected: 0.00318347509113001 ) 
##                                        0.008201701
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(pa_obj_log, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(pa_obj_log, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(pa_obj_log, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(pa_obj_log, variable = priorsense_vars)
}
# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_obj <- summarize_brms(
  pa_obj_log, 
  stats_to_report = stats_to_report,
  rope_range = rope_range_log,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) 
## Sampling priors, please wait...
## Warning in summarize_brms(pa_obj_log, stats_to_report = stats_to_report, :
## Coefficients were exponentiated. Double check if this was intended.
summary_pa_obj %>%
  print_df(rows_to_pack = rows_to_pack)
exp(Est.) SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 110.43*** 5.92 [99.31, 123.00] 1.000 [0.94, 1.07] 0.000 >100 Overwhelming Evidence 1.000 7621 15326
Within-Person Effects
Daily individual’s experienced persuasion 1.03 0.02 [ 1.00, 1.06] 0.955 [0.94, 1.07] 0.995 0.024 Very Strong Evidence for Null 1.000 30820 29811
Daily partner’s experienced persuasion 1.02 0.02 [ 0.98, 1.05] 0.833 [0.94, 1.07] 0.998 0.009 Very Strong Evidence for Null 1.000 34244 32523
Daily individual’s experienced pressure 0.94 0.03 [ 0.88, 1.01] 0.953 [0.94, 1.07] 0.583 0.020 Very Strong Evidence for Null 1.000 57050 32240
Daily partner’s experienced pressure 0.98 0.03 [ 0.92, 1.05] 0.704 [0.94, 1.07] 0.908 0.005 Very Strong Evidence for Null 1.000 66267 31179
Daily individual’s experienced pushing 1.02 0.03 [ 0.96, 1.07] 0.744 [0.94, 1.07] 0.969 0.007 Very Strong Evidence for Null 1.000 39949 31062
Daily partner’s experienced pushing 1.00 0.02 [ 0.96, 1.05] 0.544 [0.94, 1.07] 0.996 0.004 Very Strong Evidence for Null 1.000 62526 31324
Day 0.97 0.03 [ 0.91, 1.04] 0.766 [0.94, 1.07] 0.859 0.005 Very Strong Evidence for Null 1.000 85145 28779
Own Actionplan 1.07** 0.03 [ 1.02, 1.12] 0.996 [0.94, 1.07] 0.473 0.139 Moderate Evidence for Null 1.000 85924 29740
Partner Actionplan 1.05 0.03 [ 1.00, 1.10] 0.971 [0.94, 1.07] 0.753 0.024 Very Strong Evidence for Null 1.000 80280 28316
Daily weartime 1.00*** 0.00 [ 1.00, 1.00] 1.000 [0.94, 1.07] 1.000 16.395 Strong Evidence 1.000 89950 27229
Between-Person Effects
Mean individual’s experienced persuasion 1.10 0.15 [ 0.83, 1.47] 0.760 [0.94, 1.07] 0.284 0.017 Very Strong Evidence for Null 1.000 8046 14442
Mean partner’s experienced persuasion 0.99 0.14 [ 0.74, 1.32] 0.535 [0.94, 1.07] 0.350 0.013 Very Strong Evidence for Null 1.000 7934 14021
Mean individual’s experienced pressure 0.97 0.14 [ 0.73, 1.30] 0.577 [0.94, 1.07] 0.337 0.008 Very Strong Evidence for Null 1.000 10912 20180
Mean partner’s experienced pressure 0.97 0.14 [ 0.73, 1.28] 0.595 [0.94, 1.07] 0.343 0.009 Very Strong Evidence for Null 1.000 9972 18919
Mean individual’s experienced pushing 0.96 0.20 [ 0.64, 1.45] 0.587 [0.94, 1.07] 0.244 0.011 Very Strong Evidence for Null 1.000 9512 15959
Mean partner’s experienced pushing 1.21 0.24 [ 0.81, 1.83] 0.832 [0.94, 1.07] 0.159 0.018 Very Strong Evidence for Null 1.000 9390 15940
Mean weartime 1.00 0.00 [ 1.00, 1.00] 0.924 [0.94, 1.07] 1.000 0.019 Very Strong Evidence for Null 1.000 43367 31436
Random Effects
sd(Intercept) 0.30 0.04 [0.23, 0.39] NA NA NA NA NA 1.001 9175 18817
sd(Daily individual’s experienced persuasion) 0.05 0.01 [0.02, 0.08] NA NA NA NA NA 1.000 18218 14283
sd(Daily partner’s experienced persuasion) 0.06 0.02 [0.03, 0.09] NA NA NA NA NA 1.000 19503 16083
sd(Daily individual’s experienced pressure) 0.04 0.03 [0.00, 0.13] NA NA NA NA NA 1.000 21550 21109
sd(Daily partner’s experienced pressure) 0.03 0.03 [0.00, 0.11] NA NA NA NA NA 1.000 26998 21329
sd(Daily individual’s experienced pushing) 0.08 0.04 [0.01, 0.16] NA NA NA NA NA 1.000 9406 9602
sd(Daily partner’s experienced pushing) 0.03 0.03 [0.00, 0.09] NA NA NA NA NA 1.000 16849 20307
Additional Parameters
sigma 0.57 0.01 [0.56, 0.59] NA NA NA NA NA 1.000 72383 28323 
plot(
  bayestestR::p_direction(pa_obj_log),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(pa_obj_log, range = rope_range_log, ci = 1)
) + theme_bw()
## Possible multicollinearity between b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.89), b_pressure_self_cb and
##   b_persuasion_self_cb (r = 0.71), b_pressure_partner_cb and
##   b_persuasion_self_cb (r = 0.72), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.75), b_pushing_partner_cb and
##   b_pushing_self_cb (r = 0.83). This might lead to inappropriate results.
##   See 'Details' in '?rope'.

# Nothing significant, no plots

Affect

range(df_double$aff, na.rm = T)
## [1] 0 5
hist(df_double$aff, breaks = 15)

Gaussian

formula <- bf(
  aff ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    plan_self + plan_partner +
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b")
  ,brms::set_prior("normal(0, 20)", class = "Intercept", lb=1, ub=6)
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = gaussian()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

mood_gauss <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("mood_gauss", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
mood_gauss_digest <- digest::digest(mood_gauss)
if (check_models) {
  check_brms(mood_gauss, log_pp_check = FALSE)
  DHARMa.check_brms.all(mood_gauss, integer = FALSE)
}
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                   Term  VIF     VIF 95% CI Increased SE Tolerance
##     persuasion_self_cw 1.13 [ 1.10,  1.17]         1.06      0.89
##  persuasion_partner_cw 1.09 [ 1.06,  1.13]         1.04      0.92
##       pressure_self_cw 1.09 [ 1.06,  1.13]         1.04      0.92
##    pressure_partner_cw 1.07 [ 1.05,  1.12]         1.04      0.93
##        pushing_self_cw 1.17 [ 1.13,  1.21]         1.08      0.86
##     pushing_partner_cw 1.14 [ 1.11,  1.19]         1.07      0.87
##              plan_self 1.29 [ 1.25,  1.34]         1.14      0.78
##           plan_partner 1.28 [ 1.24,  1.33]         1.13      0.78
##                    day 1.05 [ 1.02,  1.09]         1.02      0.95
##  Tolerance 95% CI
##      [0.85, 0.91]
##      [0.88, 0.94]
##      [0.89, 0.95]
##      [0.90, 0.96]
##      [0.83, 0.88]
##      [0.84, 0.90]
##      [0.75, 0.80]
##      [0.75, 0.81]
##      [0.91, 0.98]
## 
## Moderate Correlation
## 
##                 Term  VIF     VIF 95% CI Increased SE Tolerance
##     pressure_self_cb 5.59 [ 5.29,  5.90]         2.36      0.18
##  pressure_partner_cb 5.46 [ 5.17,  5.76]         2.34      0.18
##      pushing_self_cb 8.63 [ 8.16,  9.12]         2.94      0.12
##   pushing_partner_cb 8.66 [ 8.20,  9.16]         2.94      0.12
##  Tolerance 95% CI
##      [0.17, 0.19]
##      [0.17, 0.19]
##      [0.11, 0.12]
##      [0.11, 0.12]
## 
## High Correlation
## 
##                   Term   VIF     VIF 95% CI Increased SE Tolerance
##     persuasion_self_cb 13.82 [13.05, 14.63]         3.72      0.07
##  persuasion_partner_cb 13.75 [12.99, 14.56]         3.71      0.07
##  Tolerance 95% CI
##      [0.07, 0.08]
##      [0.07, 0.08]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        cauchy         44%
##        normal         41%
##   exponential          6%
## 
## Predicted Distribution of Response
## 
##   Distribution Probability
##        tweedie         41%
##  beta-binomial         38%
##    half-cauchy          6%
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  model.check
## outliers at both margin(s) = 30, observations = 3684, p-value =
## 2.98e-10
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.005500816 0.011604846
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                          0.008143322
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(mood_gauss, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(mood_gauss, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(mood_gauss, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(mood_gauss, variable = priorsense_vars)
}
summary_mood <- summarize_brms(
  mood_gauss, 
  stats_to_report = stats_to_report,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = F) 
## Warning in compute_rope(range = rope_range): Collinearity detected. Some VIFs
## are > 10. This may invalidate ROPE inferences!

Check for Multicollinearity

Low Correlation

              Term  VIF     VIF 95% CI Increased SE Tolerance
persuasion_self_cw 1.13 [ 1.10,  1.17]         1.06      0.89

persuasion_partner_cw 1.09 [ 1.06, 1.13] 1.04 0.92 pressure_self_cw 1.09 [ 1.06, 1.13] 1.04 0.92 pressure_partner_cw 1.07 [ 1.05, 1.12] 1.04 0.93 pushing_self_cw 1.17 [ 1.13, 1.21] 1.08 0.86 pushing_partner_cw 1.14 [ 1.11, 1.19] 1.07 0.87 plan_self 1.29 [ 1.25, 1.34] 1.14 0.78 plan_partner 1.28 [ 1.24, 1.33] 1.13 0.78 day 1.05 [ 1.02, 1.09] 1.02 0.95 Tolerance 95% CI [0.85, 0.91] [0.88, 0.94] [0.89, 0.95] [0.90, 0.96] [0.83, 0.88] [0.84, 0.90] [0.75, 0.80] [0.75, 0.81] [0.91, 0.98]

Moderate Correlation

            Term  VIF     VIF 95% CI Increased SE Tolerance
pressure_self_cb 5.59 [ 5.29,  5.90]         2.36      0.18

pressure_partner_cb 5.46 [ 5.17, 5.76] 2.34 0.18 pushing_self_cb 8.63 [ 8.16, 9.12] 2.94 0.12 pushing_partner_cb 8.66 [ 8.20, 9.16] 2.94 0.12 Tolerance 95% CI [0.17, 0.19] [0.17, 0.19] [0.11, 0.12] [0.11, 0.12]

High Correlation

              Term   VIF     VIF 95% CI Increased SE Tolerance
persuasion_self_cb 13.82 [13.05, 14.63]         3.72      0.07

persuasion_partner_cb 13.75 [12.99, 14.56] 3.71 0.07 Tolerance 95% CI [0.07, 0.08] [0.07, 0.08]

## Sampling priors, please wait...
summary_mood %>%
  print_df(rows_to_pack = rows_to_pack)
Est. SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 3.67*** 0.10 [ 3.45, 3.87] 1.000 [-0.11, 0.11] 0.000 >100 Overwhelming Evidence 1.001 4587 9639
Within-Person Effects
Daily individual’s experienced persuasion 0.00 0.02 [-0.04, 0.04] 0.513 [-0.11, 0.11] 1.000 0.004 Very Strong Evidence for Null 1.000 41022 28769
Daily partner’s experienced persuasion 0.02 0.02 [-0.03, 0.07] 0.798 [-0.11, 0.11] 1.000 0.006 Very Strong Evidence for Null 1.000 33612 29443
Daily individual’s experienced pressure -0.02 0.05 [-0.14, 0.08] 0.669 [-0.11, 0.11] 0.944 0.004 Very Strong Evidence for Null 1.000 42267 26844
Daily partner’s experienced pressure -0.01 0.05 [-0.13, 0.10] 0.582 [-0.11, 0.11] 0.951 0.004 Very Strong Evidence for Null 1.000 39334 27679
Daily individual’s experienced pushing 0.01 0.03 [-0.06, 0.07] 0.560 [-0.11, 0.11] 0.999 0.004 Very Strong Evidence for Null 1.000 44109 30412
Daily partner’s experienced pushing 0.07* 0.03 [ 0.00, 0.14] 0.976 [-0.11, 0.11] 0.883 0.033 Strong Evidence for Null 1.000 29114 28247
Day 0.25*** 0.06 [ 0.14, 0.36] 1.000 [-0.11, 0.11] 0.006 27.512 Strong Evidence 1.000 63424 28476
Own Actionplan 0.11** 0.04 [ 0.04, 0.19] 0.998 [-0.11, 0.11] 0.502 0.243 Moderate Evidence for Null 1.000 71905 28802
Partner Actionplan -0.04 0.04 [-0.11, 0.04] 0.822 [-0.11, 0.11] 0.976 0.005 Very Strong Evidence for Null 1.000 62756 30927
Daily weartime NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion 0.34 0.28 [-0.20, 0.90] 0.894 [-0.11, 0.11] 0.154 0.028 Very Strong Evidence for Null 1.000 5550 10810
Mean partner’s experienced persuasion 0.21 0.27 [-0.33, 0.77] 0.785 [-0.11, 0.11] 0.239 0.018 Very Strong Evidence for Null 1.000 5463 10623
Mean individual’s experienced pressure -0.30 0.27 [-0.83, 0.24] 0.861 [-0.11, 0.11] 0.187 0.015 Very Strong Evidence for Null 1.000 6359 13148
Mean partner’s experienced pressure -0.32 0.27 [-0.85, 0.22] 0.877 [-0.11, 0.11] 0.174 0.016 Very Strong Evidence for Null 1.000 6160 12060
Mean individual’s experienced pushing 0.19 0.38 [-0.59, 0.95] 0.687 [-0.11, 0.11] 0.207 0.012 Very Strong Evidence for Null 1.001 6110 12468
Mean partner’s experienced pushing 0.35 0.39 [-0.43, 1.11] 0.817 [-0.11, 0.11] 0.154 0.017 Very Strong Evidence for Null 1.001 6082 10909
Mean weartime NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.60 0.07 [0.48, 0.78] NA NA NA NA NA 1.001 7253 13879
sd(Daily individual’s experienced persuasion) 0.04 0.03 [0.00, 0.10] NA NA NA NA NA 1.000 11219 15489
sd(Daily partner’s experienced persuasion) 0.07 0.03 [0.01, 0.13] NA NA NA NA NA 1.000 10416 8098
sd(Daily individual’s experienced pressure) 0.07 0.06 [0.00, 0.24] NA NA NA NA NA 1.000 16225 18606
sd(Daily partner’s experienced pressure) 0.08 0.07 [0.00, 0.26] NA NA NA NA NA 1.000 16048 19199
sd(Daily individual’s experienced pushing) 0.06 0.04 [0.00, 0.15] NA NA NA NA NA 1.000 14539 16036
sd(Daily partner’s experienced pushing) 0.08 0.05 [0.01, 0.17] NA NA NA NA NA 1.000 12571 11771
Additional Parameters
sigma 0.96 0.01 [0.94, 0.98] NA NA NA NA NA 1.000 62396 29720 
plot(
  bayestestR::p_direction(mood_gauss),
  priors = TRUE
)  + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(mood_gauss, ci = 1)
) + theme_bw()
## Possible multicollinearity between b_pressure_self_cb and
##   b_persuasion_self_cb (r = 0.8), b_pressure_partner_cb and
##   b_persuasion_self_cb (r = 0.77), b_pressure_self_cb and
##   b_persuasion_partner_cb (r = 0.78), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.8), b_pressure_partner_cb and
##   b_pressure_self_cb (r = 0.75), b_pushing_partner_cb and
##   b_pushing_self_cb (r = 0.89). This might lead to inappropriate results.
##   See 'Details' in '?rope'.

conditional_spaghetti(
  mood_gauss, 
  effects = c('pushing_partner_cw'),
  group_var = 'coupleID',
  plot_full_range = TRUE
)

$pushing_partner_cw

Reactance

range(df_double$reactance, na.rm = T) 
## [1] 0 5
hist(df_double$reactance, breaks = 7) 

hist(log(df_double$reactance+0.1), breaks = 10)

Ordinal

df_double$reactance_ordinal <- factor(df_double$reactance,
                                      levels = 0:5, 
                                      ordered = TRUE)

formula <- bf(
  reactance_ordinal ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    plan_self + plan_partner +
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = cumulative() # HURDLE_CUMULATIVE
#  )


#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

reactance_ordinal <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::cumulative(),
  control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777
  , file = file.path("models_cache_brms", paste0("reactance_ordinal", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
reactance_ordinal_digest <- digest::digest(reactance_ordinal)
if (check_models) {
  check_brms(reactance_ordinal)
  DHARMa.check_brms.all(reactance_ordinal, outliers_type = 'bootstrap')
}
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                   Term  VIF    VIF 95% CI Increased SE Tolerance
##       pressure_self_cw 4.44 [4.13,  4.78]         2.11      0.23
##    pressure_partner_cw 1.50 [1.43,  1.59]         1.23      0.66
##        pushing_self_cw 1.30 [1.24,  1.37]         1.14      0.77
##     pushing_partner_cw 1.12 [1.07,  1.18]         1.06      0.90
##     persuasion_self_cb 1.09 [1.05,  1.15]         1.04      0.92
##  persuasion_partner_cb 1.08 [1.04,  1.14]         1.04      0.93
##       pressure_self_cb 1.39 [1.32,  1.47]         1.18      0.72
##    pressure_partner_cb 1.19 [1.14,  1.26]         1.09      0.84
##                    day 2.13 [2.00,  2.27]         1.46      0.47
##  Tolerance 95% CI
##      [0.21, 0.24]
##      [0.63, 0.70]
##      [0.73, 0.81]
##      [0.85, 0.93]
##      [0.87, 0.95]
##      [0.87, 0.96]
##      [0.68, 0.76]
##      [0.79, 0.88]
##      [0.44, 0.50]
## 
## Moderate Correlation
## 
##                Term  VIF    VIF 95% CI Increased SE Tolerance Tolerance 95% CI
##  persuasion_self_cw 8.92 [8.25,  9.64]         2.99      0.11     [0.10, 0.12]
##     pushing_self_cb 5.87 [5.45,  6.33]         2.42      0.17     [0.16, 0.18]
##  pushing_partner_cb 7.14 [6.61,  7.71]         2.67      0.14     [0.13, 0.15]
##           plan_self 5.60 [5.19,  6.04]         2.37      0.18     [0.17, 0.19]
##        plan_partner 6.19 [5.74,  6.68]         2.49      0.16     [0.15, 0.17]
## 
## High Correlation
## 
##                   Term   VIF    VIF 95% CI Increased SE Tolerance
##  persuasion_partner_cw 10.29 [9.51, 11.13]         3.21      0.10
##  Tolerance 95% CI
##      [0.09, 0.11]
## Error in h(simpleError(msg, call)) : 
##   error in evaluating the argument 'x' in selecting a method for function 'print': Predictive errors are not defined for ordinal or categorical models.
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 3 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 0, observations = 737, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.002713704
## sample estimates:
## outlier frequency (expected: 0.00039348710990502 ) 
##                                                  0
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(reactance_ordinal, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(reactance_ordinal, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(reactance_ordinal, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(reactance_ordinal, variable = priorsense_vars)
}
summary_reactance_ordinal <- summarize_brms(
  reactance_ordinal, 
  stats_to_report = stats_to_report,
  rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed_ordinal,
  model_rows_random = model_rows_random_ordinal,
  model_rownames_fixed = model_rownames_fixed_ordinal,
  model_rownames_random = model_rownames_random_ordinal,
  exponentiate = T) 
## Warning in compute_rope(range = rope_range): Collinearity detected. Some VIFs
## are > 10. This may invalidate ROPE inferences!

Check for Multicollinearity

Low Correlation

              Term  VIF    VIF 95% CI Increased SE Tolerance
  pressure_self_cw 4.44 [4.13,  4.78]         2.11      0.23

pressure_partner_cw 1.50 [1.43, 1.59] 1.23 0.66 pushing_self_cw 1.30 [1.24, 1.37] 1.14 0.77 pushing_partner_cw 1.12 [1.07, 1.18] 1.06 0.90 persuasion_self_cb 1.09 [1.05, 1.15] 1.04 0.92 persuasion_partner_cb 1.08 [1.04, 1.14] 1.04 0.93 pressure_self_cb 1.39 [1.32, 1.47] 1.18 0.72 pressure_partner_cb 1.19 [1.14, 1.26] 1.09 0.84 day 2.13 [2.00, 2.27] 1.46 0.47 Tolerance 95% CI [0.21, 0.24] [0.63, 0.70] [0.73, 0.81] [0.85, 0.93] [0.87, 0.95] [0.87, 0.96] [0.68, 0.76] [0.79, 0.88] [0.44, 0.50]

Moderate Correlation

           Term  VIF    VIF 95% CI Increased SE Tolerance Tolerance 95% CI

persuasion_self_cw 8.92 [8.25, 9.64] 2.99 0.11 [0.10, 0.12] pushing_self_cb 5.87 [5.45, 6.33] 2.42 0.17 [0.16, 0.18] pushing_partner_cb 7.14 [6.61, 7.71] 2.67 0.14 [0.13, 0.15] plan_self 5.60 [5.19, 6.04] 2.37 0.18 [0.17, 0.19] plan_partner 6.19 [5.74, 6.68] 2.49 0.16 [0.15, 0.17]

High Correlation

              Term   VIF    VIF 95% CI Increased SE Tolerance

persuasion_partner_cw 10.29 [9.51, 11.13] 3.21 0.10 Tolerance 95% CI [0.09, 0.11]

## Sampling priors, please wait...
summary_reactance_ordinal %>%
  print_df(rows_to_pack = rows_to_pack_ordinal)
OR SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept NA NA NA NA NA NA NA NA NA NA NA
Intercept[1] 3.26*** 1.00 [ 1.81, 6.11] 1.000 [0.84, 1.20] 0.001 >100 Overwhelming Evidence 1.000 49874 34935
Intercept[2] 7.14*** 2.26 [ 3.90, 13.65] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 50185 33484
Intercept[3] 19.20*** 6.54 [ 10.08, 38.31] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 52162 35394
Intercept[4] 81.68*** 31.62 [ 39.02, 180.37] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 56054 33021
Intercept[5] 2436.87*** 1607.60 [735.28, 10016.43] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 72041 30135
Within-Person Effects
Daily individual’s experienced persuasion 0.85* 0.07 [ 0.71, 1.00] 0.978 [0.84, 1.20] 0.559 0.384 Weak Evidence for Null 1.000 51294 33052
Daily partner’s experienced persuasion 1.02 0.10 [ 0.83, 1.24] 0.566 [0.84, 1.20] 0.920 0.045 Strong Evidence for Null 1.000 48711 32207
Daily individual’s experienced pressure 1.84* 0.35 [ 1.17, 2.64] 0.993 [0.84, 1.20] 0.029 2.297 Weak Evidence 1.000 34038 28803
Daily partner’s experienced pressure 1.23 0.30 [ 0.68, 2.15] 0.792 [0.84, 1.20] 0.384 0.081 Strong Evidence for Null 1.000 33901 24427
Daily individual’s experienced pushing 1.20 0.13 [ 0.97, 1.51] 0.951 [0.84, 1.20] 0.489 0.184 Moderate Evidence for Null 1.000 49980 34581
Daily partner’s experienced pushing 0.93 0.12 [ 0.70, 1.22] 0.700 [0.84, 1.20] 0.759 0.052 Strong Evidence for Null 1.000 47206 30314
Day 1.56 0.56 [ 0.78, 3.10] 0.891 [0.84, 1.20] 0.193 0.081 Strong Evidence for Null 1.000 74170 33770
Own Actionplan 0.83 0.24 [ 0.47, 1.47] 0.743 [0.84, 1.20] 0.384 0.049 Strong Evidence for Null 1.000 78721 30650
Partner Actionplan 0.92 0.24 [ 0.55, 1.53] 0.624 [0.84, 1.20] 0.493 0.040 Strong Evidence for Null 1.000 75874 29696
Daily weartime NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion 0.99 0.53 [ 0.33, 2.98] 0.505 [0.84, 1.20] 0.263 0.058 Strong Evidence for Null 1.000 33318 32341
Mean partner’s experienced persuasion 1.46 0.90 [ 0.44, 5.10] 0.732 [0.84, 1.20] 0.195 0.067 Strong Evidence for Null 1.000 36161 30998
Mean individual’s experienced pressure 3.86* 2.23 [ 1.25, 12.59] 0.990 [0.84, 1.20] 0.017 0.794 Weak Evidence for Null 1.000 30585 33799
Mean partner’s experienced pressure 1.21 0.73 [ 0.35, 3.88] 0.622 [0.84, 1.20] 0.223 0.054 Strong Evidence for Null 1.000 31105 31607
Mean individual’s experienced pushing 1.16 0.91 [ 0.25, 5.75] 0.575 [0.84, 1.20] 0.177 0.057 Strong Evidence for Null 1.000 28869 30895
Mean partner’s experienced pushing 0.09* 0.09 [ 0.01, 0.57] 0.994 [0.84, 1.20] 0.006 1.572 Weak Evidence 1.000 35396 32561
Mean weartime NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.80 0.20 [0.46, 1.28] NA NA NA NA NA 1.000 19098 28104
sd(Daily individual’s experienced persuasion) 0.18 0.13 [0.01, 0.44] NA NA NA NA NA 1.001 9634 16299
sd(Daily partner’s experienced persuasion) 0.21 0.14 [0.01, 0.53] NA NA NA NA NA 1.000 13095 16639
sd(Daily individual’s experienced pressure) 0.53 0.26 [0.07, 1.14] NA NA NA NA NA 1.000 11427 12889
sd(Daily partner’s experienced pressure) 0.44 0.41 [0.02, 1.59] NA NA NA NA NA 1.000 12940 23080
sd(Daily individual’s experienced pushing) 0.21 0.14 [0.01, 0.52] NA NA NA NA NA 1.000 14083 18971
sd(Daily partner’s experienced pushing) 0.16 0.15 [0.01, 0.63] NA NA NA NA NA 1.000 20365 25999
Additional Parameters
sigma NA NA NA NA NA NA NA NA NA NA NA
disc 1.00 0.00 [1.00, 1.00] NA NA NA NA NA NA NA NA 
plot(
  bayestestR::p_direction(reactance_ordinal),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(reactance_ordinal, range = c(-0.18, 0.18), ci = 1)
) + theme_bw()
## Possible multicollinearity between b_Intercept[4] and b_Intercept[2] (r
##   = 0.81), b_Intercept[4] and b_Intercept[3] (r = 0.86),
##   b_pressure_self_cb and b_persuasion_self_cb (r = 0.71),
##   b_pressure_partner_cb and b_persuasion_partner_cb (r = 0.79). This might
##   lead to inappropriate results. See 'Details' in '?rope'.

conditional_spaghetti(
  reactance_ordinal, 
  effects = c('persuasion_self_cw', 'pressure_self_cw')
  , group_var = 'coupleID'
  #, n_groups = 15
  , plot_full_range = T
)

\(persuasion_self_cw <img src="01_FinalModelsPlan_files/figure-html/report_reactance_ordinal-3.png" width="2400" />\)pressure_self_cw

Binary

introduce_binary_reactance <- function(data) {
  data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
  return(data)
}



df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
  for (i in seq_along(implist)) {
    implist[[i]] <- introduce_binary_reactance(implist[[i]])
  }
}


formula <- bf(
  is_reactance ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    plan_self + plan_partner +
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
  )



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5) 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = bernoulli()
#  )



#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

is_reactance <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::bernoulli(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("is_reactance", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
is_reactance_digest <- digest::digest(is_reactance)
if (check_models) {
  check_brms(is_reactance)
  DHARMa.check_brms.all(is_reactance, integer = FALSE)
}
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                   Term  VIF   VIF 95% CI Increased SE Tolerance
##     persuasion_self_cw 1.06 [1.03, 1.13]         1.03      0.94
##  persuasion_partner_cw 1.17 [1.12, 1.23]         1.08      0.86
##       pressure_self_cw 1.05 [1.02, 1.12]         1.02      0.96
##    pressure_partner_cw 1.05 [1.02, 1.13]         1.03      0.95
##        pushing_self_cw 1.25 [1.19, 1.32]         1.12      0.80
##     pushing_partner_cw 1.16 [1.11, 1.22]         1.08      0.86
##     persuasion_self_cb 2.78 [2.60, 2.98]         1.67      0.36
##  persuasion_partner_cb 3.34 [3.12, 3.59]         1.83      0.30
##       pressure_self_cb 3.00 [2.80, 3.22]         1.73      0.33
##    pressure_partner_cb 3.07 [2.86, 3.29]         1.75      0.33
##        pushing_self_cb 2.27 [2.13, 2.43]         1.51      0.44
##     pushing_partner_cb 2.11 [1.99, 2.26]         1.45      0.47
##              plan_self 1.67 [1.58, 1.77]         1.29      0.60
##           plan_partner 1.55 [1.47, 1.65]         1.25      0.64
##                    day 1.06 [1.03, 1.13]         1.03      0.94
##  Tolerance 95% CI
##      [0.88, 0.97]
##      [0.81, 0.89]
##      [0.89, 0.98]
##      [0.89, 0.98]
##      [0.76, 0.84]
##      [0.82, 0.90]
##      [0.34, 0.38]
##      [0.28, 0.32]
##      [0.31, 0.36]
##      [0.30, 0.35]
##      [0.41, 0.47]
##      [0.44, 0.50]
##      [0.56, 0.63]
##      [0.61, 0.68]
##      [0.88, 0.97]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##       tweedie         34%
##        normal         19%
##     bernoulli          9%
## 
## Predicted Distribution of Response
## 
##   Distribution Probability
##      bernoulli         78%
##  beta-binomial         12%
##       binomial          9%
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 40 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  model.check
## outliers at both margin(s) = 0, observations = 737, p-value = 0.4132
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.000000000 0.004992758
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                                    0
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(is_reactance, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(is_reactance, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(is_reactance, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(is_reactance, variable = priorsense_vars)
}
summary_is_reactance <- summarize_brms(
  is_reactance, 
  stats_to_report = stats_to_report,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) 
## Sampling priors, please wait...
summary_is_reactance %>%
  print_df(rows_to_pack = rows_to_pack)
OR SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 0.37** 0.12 [0.19, 0.72] 0.998 [0.83, 1.20] 0.008 3.304 Moderate Evidence 1.000 31019 30682
Within-Person Effects
Daily individual’s experienced persuasion 0.84 0.08 [0.68, 1.02] 0.962 [0.83, 1.20] 0.517 0.293 Moderate Evidence for Null 1.000 32119 29620
Daily partner’s experienced persuasion 1.12 0.17 [0.83, 1.57] 0.776 [0.83, 1.20] 0.648 0.090 Strong Evidence for Null 1.000 22128 24649
Daily individual’s experienced pressure 2.02* 0.67 [1.04, 4.77] 0.979 [0.83, 1.20] 0.047 1.064 Weak Evidence 1.001 19360 18037
Daily partner’s experienced pressure 1.47 0.68 [0.55, 4.90] 0.804 [0.83, 1.20] 0.218 0.142 Moderate Evidence for Null 1.000 17715 17902
Daily individual’s experienced pushing 1.33* 0.18 [1.03, 1.77] 0.985 [0.83, 1.20] 0.215 0.592 Weak Evidence for Null 1.000 31430 29380
Daily partner’s experienced pushing 0.89 0.18 [0.58, 1.35] 0.719 [0.83, 1.20] 0.557 0.083 Strong Evidence for Null 1.000 29538 26005
Day 1.71 0.70 [0.77, 3.81] 0.907 [0.83, 1.20] 0.151 0.108 Moderate Evidence for Null 1.000 41465 30365
Own Actionplan 0.86 0.28 [0.45, 1.63] 0.681 [0.83, 1.20] 0.382 0.050 Strong Evidence for Null 1.000 42283 30081
Partner Actionplan 0.86 0.26 [0.48, 1.56] 0.686 [0.83, 1.20] 0.407 0.049 Strong Evidence for Null 1.000 41935 29405
Daily weartime NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion 1.82 1.26 [0.47, 7.74] 0.810 [0.83, 1.20] 0.144 0.108 Moderate Evidence for Null 1.000 17504 23517
Mean partner’s experienced persuasion 2.08 1.60 [0.48, 10.26] 0.834 [0.83, 1.20] 0.125 0.109 Moderate Evidence for Null 1.000 19179 25934
Mean individual’s experienced pressure 34.42** 44.66 [3.20, 583.15] 0.999 [0.83, 1.20] 0.002 10.092 Strong Evidence 1.000 13676 20206
Mean partner’s experienced pressure 1.64 2.21 [0.09, 20.83] 0.644 [0.83, 1.20] 0.102 0.122 Moderate Evidence for Null 1.000 13260 22817
Mean individual’s experienced pushing 0.56 0.65 [0.05, 5.49] 0.692 [0.83, 1.20] 0.110 0.093 Strong Evidence for Null 1.000 15142 22421
Mean partner’s experienced pushing 0.06* 0.07 [0.00, 0.67] 0.988 [0.83, 1.20] 0.008 1.070 Weak Evidence 1.000 15687 23090
Mean weartime NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 1.22 0.27 [0.78, 1.86] NA NA NA NA NA 1.000 13237 20425
sd(Daily individual’s experienced persuasion) 0.24 0.15 [0.02, 0.55] NA NA NA NA NA 1.000 7747 11825
sd(Daily partner’s experienced persuasion) 0.51 0.21 [0.12, 1.03] NA NA NA NA NA 1.000 10640 9083
sd(Daily individual’s experienced pressure) 1.08 0.57 [0.12, 2.46] NA NA NA NA NA 1.002 6773 7302
sd(Daily partner’s experienced pressure) 0.91 0.73 [0.05, 2.91] NA NA NA NA NA 1.000 11490 15623
sd(Daily individual’s experienced pushing) 0.26 0.17 [0.02, 0.64] NA NA NA NA NA 1.000 11614 13089
sd(Daily partner’s experienced pushing) 0.28 0.25 [0.01, 1.03] NA NA NA NA NA 1.001 13159 17480
Additional Parameters
sigma NA NA NA NA NA NA NA NA NA NA NA 
plot(
  bayestestR::p_direction(is_reactance),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(is_reactance, ci = 1)
) + theme_bw()

conditional_spaghetti(
  is_reactance, 
  effects = c('pressure_self_cw', 'pushing_self_cw'),
  group_var = 'coupleID',
  plot_full_range = TRUE
)

\(pressure_self_cw <img src="01_FinalModelsPlan_files/figure-html/report_is_reactance-3.png" width="2400" />\)pushing_self_cw 

hypothesis(is_reactance, "exp(pressure_self_cw) > exp(pushing_self_cw)")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (exp(pressure_sel... > 0     0.89      1.02    -0.28      2.7       6.91
##   Post.Prob Star
## 1      0.87     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Report All Models

summary_all_models <- report_side_by_side(
  pa_sub,
  pa_obj_log,
  mood_gauss,
  is_reactance,
  
  stats_to_report = c('CI'),
  model_rows_random = model_rows_random,
  model_rows_fixed = model_rows_fixed,
  model_rownames_random = model_rownames_random,
  model_rownames_fixed = model_rownames_fixed
)

[1] “pa_sub”

## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “pa_obj_log”

## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “mood_gauss” [1] “is_reactance”

summary_all_models <- summary_all_models %>%
  print_df(rows_to_pack = rows_to_pack) %>%
  add_header_above(
    c(
      " ", "Hurdle Component" = 2, "Non-Zero Component" = 2,
      " " = 6
    )
  ) %>%
  add_header_above(
    c(" ", "Subjective MVPA Hurdle Lognormal" = 4,  
      "Device-Based MVPA Log (Gaussian)" = 2, 
      "Mood Gaussian" = 2,
      #"Reactance Ordinal" = 2,
      "Reactance Dichotome" = 2
    )
  )

export_xlsx(
  summary_all_models, 
  rows_to_pack = rows_to_pack,
  file.path("Output", paste0("AllModels", suffix, ".xlsx")), 
  merge_option = 'both', 
  simplify_2nd_row = FALSE,
  line_above_rows = c(1,2),
  line_below_rows = c(-1)
)
## 
## Attaching package: 'rvest'
## The following object is masked from 'package:readr':
## 
##     guess_encoding
summary_all_models
Subjective MVPA Hurdle Lognormal
Device-Based MVPA Log (Gaussian)
Mood Gaussian
Reactance Dichotome
Hurdle Component
Non-Zero Component
exp(Est.)_hu pa_sub 95% CI_hu pa_sub exp(Est.)_nonzero pa_sub 95% CI_nonzero pa_sub exp(Est.) pa_obj_log 95% CI pa_obj_log Est. mood_gauss 95% CI mood_gauss OR is_reactance 95% CI is_reactance
Intercept 0.21*** [ 0.15, 0.30] 35.17*** [29.76, 41.52] 110.43*** [99.31, 123.00] 3.67*** [ 3.45, 3.87] 0.37** [0.19, 0.72]
Within-Person Effects
Daily individual’s experienced persuasion 1.58*** [ 1.37, 1.84] 1.03 [ 0.98, 1.09] 1.03 [ 1.00, 1.06] 0.00 [-0.04, 0.04] 0.84 [0.68, 1.02]
Daily partner’s experienced persuasion 1.34*** [ 1.18, 1.53] 1.03 [ 0.99, 1.08] 1.02 [ 0.98, 1.05] 0.02 [-0.03, 0.07] 1.12 [0.83, 1.57]
Daily individual’s experienced pressure 0.98 [ 0.69, 1.36] 0.90* [ 0.80, 0.99] 0.94 [ 0.88, 1.01] -0.02 [-0.14, 0.08] 2.02* [1.04, 4.77]
Daily partner’s experienced pressure 1.55* [ 1.04, 2.67] 0.95 [ 0.86, 1.03] 0.98 [ 0.92, 1.05] -0.01 [-0.13, 0.10] 1.47 [0.55, 4.90]
Daily individual’s experienced pushing 1.01 [ 0.76, 1.37] 0.99 [ 0.93, 1.06] 1.02 [ 0.96, 1.07] 0.01 [-0.06, 0.07] 1.33* [1.03, 1.77]
Daily partner’s experienced pushing 1.34** [ 1.07, 1.72] 0.96 [ 0.90, 1.02] 1.00 [ 0.96, 1.05] 0.07* [ 0.00, 0.14] 0.89 [0.58, 1.35]
Day 0.81 [ 0.60, 1.09] 0.98 [ 0.87, 1.11] 0.97 [ 0.91, 1.04] 0.25*** [ 0.14, 0.36] 1.71 [0.77, 3.81]
Own Actionplan 10.50*** [ 8.57, 12.90] 1.34*** [ 1.21, 1.47] 1.07** [ 1.02, 1.12] 0.11** [ 0.04, 0.19] 0.86 [0.45, 1.63]
Partner Actionplan 1.17 [ 0.95, 1.43] 1.07 [ 0.98, 1.17] 1.05 [ 1.00, 1.10] -0.04 [-0.11, 0.04] 0.86 [0.48, 1.56]
Daily weartime NA NA NA NA 1.00*** [ 1.00, 1.00] NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion 1.41 [ 0.65, 3.02] 0.96 [ 0.70, 1.30] 1.10 [ 0.83, 1.47] 0.34 [-0.20, 0.90] 1.82 [0.47, 7.74]
Mean partner’s experienced persuasion 1.32 [ 0.61, 2.80] 0.90 [ 0.66, 1.22] 0.99 [ 0.74, 1.32] 0.21 [-0.33, 0.77] 2.08 [0.48, 10.26]
Mean individual’s experienced pressure 0.17*** [ 0.05, 0.49] 0.76 [ 0.36, 1.63] 0.97 [ 0.73, 1.30] -0.30 [-0.83, 0.24] 34.42** [3.20, 583.15]
Mean partner’s experienced pressure 0.25** [ 0.08, 0.72] 0.56 [ 0.26, 1.20] 0.97 [ 0.73, 1.28] -0.32 [-0.85, 0.22] 1.64 [0.09, 20.83]
Mean individual’s experienced pushing 1.42 [ 0.45, 4.62] 1.51 [ 0.87, 2.62] 0.96 [ 0.64, 1.45] 0.19 [-0.59, 0.95] 0.56 [0.05, 5.49]
Mean partner’s experienced pushing 2.58 [ 0.81, 8.39] 1.75* [ 1.00, 3.04] 1.21 [ 0.81, 1.83] 0.35 [-0.43, 1.11] 0.06* [0.00, 0.67]
Mean weartime NA NA NA NA 1.00 [ 1.00, 1.00] NA NA NA NA
Random Effects
sd(Intercept) 0.81 [0.62, 1.08] 0.30 [0.23, 0.40] 0.30 [0.23, 0.39] 0.60 [0.48, 0.78] 1.22 [0.78, 1.86]
sd(Daily individual’s experienced persuasion) 0.22 [0.04, 0.41] 0.11 [0.07, 0.16] 0.05 [0.02, 0.08] 0.04 [0.00, 0.10] 0.24 [0.02, 0.55]
sd(Daily partner’s experienced persuasion) 0.15 [0.01, 0.35] 0.08 [0.04, 0.13] 0.06 [0.03, 0.09] 0.07 [0.01, 0.13] 0.51 [0.12, 1.03]
sd(Daily individual’s experienced pressure) 0.20 [0.01, 0.77] 0.07 [0.00, 0.23] 0.04 [0.00, 0.13] 0.07 [0.00, 0.24] 1.08 [0.12, 2.46]
sd(Daily partner’s experienced pressure) 0.36 [0.02, 1.21] 0.06 [0.00, 0.18] 0.03 [0.00, 0.11] 0.08 [0.00, 0.26] 0.91 [0.05, 2.91]
sd(Daily individual’s experienced pushing) 0.51 [0.20, 0.93] 0.09 [0.02, 0.17] 0.08 [0.01, 0.16] 0.06 [0.00, 0.15] 0.26 [0.02, 0.64]
sd(Daily partner’s experienced pushing) 0.21 [0.01, 0.57] 0.08 [0.01, 0.16] 0.03 [0.00, 0.09] 0.08 [0.01, 0.17] 0.28 [0.01, 1.03]
Additional Parameters
sigma NA NA 0.68 [0.65, 0.70] 0.57 [0.56, 0.59] 0.96 [0.94, 0.98] NA NA 
report::report_system()

Analyses were conducted using the R Statistical language (version 4.4.2; R Core Team, 2024) on Windows 11 x64 (build 26100)

report::report_packages()
  • beepr (version 2.0; Bååth R, 2024)
  • R.methodsS3 (version 1.8.2; Bengtsson H, 2003)
  • R.oo (version 1.27.0; Bengtsson H, 2003)
  • R.utils (version 2.12.3; Bengtsson H, 2023)
  • brms (version 2.22.0; Bürkner P, 2017)
  • posterior (version 1.6.0; Bürkner P et al., 2024)
  • extrafont (version 0.19; Chang W, 2023)
  • digest (version 0.6.37; Eddelbuettel D, 2024)
  • Rcpp (version 1.0.13.1; Eddelbuettel D et al., 2024)
  • bayesplot (version 1.11.1; Gabry J, Mahr T, 2024)
  • lubridate (version 1.9.3; Grolemund G, Wickham H, 2011)
  • DHARMa (version 0.4.7; Hartig F, 2024)
  • wbCorr (version 0.1.22; Küng P, 2023)
  • see (version 0.9.0; Lüdecke D et al., 2021)
  • tibble (version 3.2.1; Müller K, Wickham H, 2023)
  • patchwork (version 1.3.0; Pedersen T, 2024)
  • R (version 4.4.2; R Core Team, 2024)
  • openxlsx (version 4.2.7.1; Schauberger P, Walker A, 2024)
  • ggplot2 (version 3.5.1; Wickham H, 2016)
  • forcats (version 1.0.0; Wickham H, 2023)
  • stringr (version 1.5.1; Wickham H, 2023)
  • rvest (version 1.0.4; Wickham H, 2024)
  • tidyverse (version 2.0.0; Wickham H et al., 2019)
  • readxl (version 1.4.3; Wickham H, Bryan J, 2023)
  • dplyr (version 1.1.4; Wickham H et al., 2023)
  • purrr (version 1.0.2; Wickham H, Henry L, 2023)
  • readr (version 2.1.5; Wickham H et al., 2024)
  • xml2 (version 1.3.6; Wickham H et al., 2023)
  • tidyr (version 1.3.1; Wickham H et al., 2024)
  • ggridges (version 0.5.6; Wilke C, 2024)
  • knitr (version 1.49; Xie Y, 2024)
  • kableExtra (version 1.4.0; Zhu H, 2024)
report::cite_packages()